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ABSTRACT 

We present new results characterizing cosmological shocks within adaptive mesh refinement N- 
Body/hydrodynamic simulations that are used to predict non-thermal components of large-scale 
structure. This represents the first study of shocks using adaptive mesh refinement. We propose 
a modified algorithm for finding shocks from those used on unigrid simulations that reduces the shock 
frequency of low Mach number shocks by a factor of ~ 3. We then apply our new technique to a large, 
(512 Mpc/h) 3 , cosmological volume and study the shock Mach number (M) distribution as a function 
of pre-shock temperature, density, and redshift. Because of the large volume of the simulation, we 
have superb statistics that results from having thousands of galaxy clusters. We find that the Mach 
number evolution can be interpreted as a method to visualize large-scale structure formation. Shocks 
with M < 5 typically trace mergers and complex flows, while 5 < A4 < 20 and Ai > 20 generally 
follow accretion onto filaments and galaxy clusters, respectively. By applying results from nonlinear 
diffusive shock acceleration models using the first-order Fermi process, we calculate the amount of 
kinetic energy that is converted into cosmic ray protons. The acceleration of cosmic ray protons is 
large enough that in order to use galaxy clusters as cosmological probes, the dynamic response of the 
gas to the cosmic rays must be included in future numerical simulations. 

Subject headings: cosmology: theory — hydrodynamics — methods: numerical — cosmic rays 



1. INTRODUCTION 

What determines the thermal history of galaxy clus- 
ters? On large scales, it is governed by the in-fall of ma- 
terial onto dark matter halos and the conversion of grav- 
itational potential energy into thermal energy. This pro- 
cess occurs through the heating of the gas via strong ac- 
cretion shocks sur r ounding galaxy clus t ers and filaments 
(iRvu et all 120031: lMiniati_e^al.l 1200 It iPfrommer et all 
120061: iPavlidou fc Fields! I2006T ). Once inside collapsed 
structures, complex flows associated with the merging 
of subhalos continue to create moderate-strength shocks 
that allow the halos to virialize. Because of this, shocks 
encode information about structure formation and its 
thermal effects on the gas. 

Cosmological shocks affect three important realms of 
structure formation and leave feedback on the surround- 
ing medium. First, shocks thermalize the incoming 
gas, providing much of the pressure support in baryons. 
This process is the basis upon which clusters are able 
to virialize. Additionally, the thermalization of gas at 
the standing accretion shocks surrounding large-scale 
filaments produces the warm-hot interclust er medium 
(WHIM) at temperatures of lO 5 ^ - 10 7 K (jDave et al.1 
1999; iCen fc OstrikerlUg gg) . The history of the mass flux 
through these shocks describes the evolution of gas in the 



Electronic address: samuel.skillman@colorado.edu 

Center for Astrophysics and Space Astronomy, Department of 
Astrophysical & Planetary Science, University of Colorado, Boul- 
der, CO 80309 

2 Theoretical Astrophysics (T-6), Los Alamos National Labora- 
tory, Los Alamos, NM 87545 

3 Department of Physics & Astronomy, Michigan State Univer- 
sity, East Lansing, MI, 48824 

4 National Science Foundation Astronomy and Astrophysics 
Postdoctoral Fellow 

5 Center for Astrophysics and Space Sciences, University of 
California-San Diego, 9500 Oilman Drive, La Jolla, CA 92093 



WHIM phase (|Pfrommer et al.ll2008ft . 

Second, the strength of the outer accretion shocks onto 
halos, characterized by the Mach number, is dependent 
upon the mass of the gravitating object. This is be- 
cause the higher mass generates larger acceleration and 
velocity in the diffuse gas while the sound speed of the 
upstream gas is uniform for previously unshocked gas, 
corresponding to a temperature of « 10 4 K. This tem- 
perature floor is created by the reionization from stars. 
Thus, the Mach number of accretion shocks can be used 
as an independent measure of cluster mass. This could 
conceivably be a powerful new tool for cluster mass esti- 
mation if we are able to observe the accretion sh ock with 
radio observations (e.g. iGiacintucci et al.|[2008f ). 

Finally, because these shocks are collisionless fea- 
tures whose interactions in the hot plasma are me- 
diated by electromagnetic fields, it is possible for a 
portion of the thermal distribution of particles to be 
accelerated and transformed into non-thermal popula- 
tions of cosmic rays (CRs) through the process of diffu- 
sive shock acceleration (P SA; e.g. iDrurv fc Falld 119861 : 
Bla ndford fc Eichleri fl98l . This process results in a 
fraction of the kinetic energy of shocking gas being con- 
verted into both t hermal and non-the r mal c omponents 
(|Kang et all 120021 : iKang fc JonesH200l [2001 . The cos- 
mic ray electron populations are likely sources of radio 
halos and radio relics i n galaxy clusters (IPfrommer et al.l 
120081: iKim et al.lll989trGiovannini et al.lll993l ). while the 
cosmic ray protons may be sources of 7-ray emission 
through their interactions with gas protons. If a sig- 
nificant portion of the gas pressure resides in cosmic 
rays, then it will likely affect gas mass fraction esti- 
mates as well as the assumption of hydrostatic equilib- 
rium. Because of the importance of these mass estimates 
in measuring dark energy, we must include the under- 
lying physics in order to perform precision cosmology 
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(| Allen et alJl2008ft . 

To date, studies of cosmological shocks have included 
observational, theoretical, and numerical techniques. 
Observationally, the majority of the work surrounding 
cosmic shocks are related to radio reli cs, of which only 
a few have been studied in d e pth (e.g. Qjottgering et al.1 
fl997t [Clarke fc Ensslinl [20061 : lOrru et all I2007D . Using 
the spectral index of the non-thermal particle distribu- 
tion, we can infer a Mach number if t he acceleration is 
due t o first order Fermi acceleration (jGiacintucci et all 
2008). Additionally, GLAST will begin observing 7-rays 
and will likely see signatures from galaxy clus ters due to 
hadro nic cosmic ray interactions with pions (Pfrommcr 
l2008h . 

On the theoretical side, the majority of analyses 
are bas ed upon manipulating the Press-Schechter for- 
malism (jPress fc Schechterlll974HSheth fe Tormenll 19991) 
to deduce first the mass function of accreting ob- 
jects and then their in t eract ions with infalling mate- 
rial. iPavlidou fc Fields! (|2006l ) extended these analyses 
to calculate the energy and mass flux through accretion 
shocks. Furthermore, several analytical attempts have 
been made to describe merger shocks, including those by 
iFuiita fc Sarazir] (pOOlh and iGabici fe Blasl (|2003allbT l . 
However, it is quite difficult to account for the complex 
morphologies that arise during structure formation using 
purely analytical frameworks. For this reason, multiple 
numerical techniques have been developed using hydro- 
dynamical simulations. 

There have been numerical studies of shocks using 
both Eulerian "sing le-grid" codes fe.g. iRvu et al.1 120031 
iKang fc Jones! l2007f ) as w ell as smoothed particle hy- 
drodynamics ( SPH) co des (jPfrommer et al.ll2006l . I2007L 
l2008t IPfVmmr^[2008h . There are advantages and dis- 
advantages of both methods. In previous work using 
grid-based codes, shocks were analyzed during post- 
processing by e xamining temper ature jumps for a given 
point in time (|Rvu et alj [2003) . However, it was im- 
possible to cover the spatial dynamic range needed to 
describe both the complex flow within halos and their 
coupling to large scale structures because of the use of 
a uniform grid. Therefore, even the largest simulations, 
with 1024 3 cells in a IQQh^Mpc volum e, has a comovin g 
spatial resolution of only 97.7h~ 1 kpc (|Rvu et al.ll2003f ). 
There have been recent attempts at prescribing hybrid 
models to study turbulent generation, but full resolution 
convergence of th e results have still not yet been achieved 
(|Rvu et al.ll2008[ ) . The advantage of a grid code is its su- 
perb shock-capturing algorithms that do not rely on the 
use of artificial viscosity when using higher-order meth- 
ods (jO'Sheaet al.ll2THT5h . 

In contrast, SPH codes are implicitly adaptive in space 
due to the Lagrangian nature of the method, e.g. high 
density regions are resolved by a larger number of par- 
ticles than low-density regions. This approach conserves 
hydrodynamic quantities exactly when they are advected 
with the flow. However, because the properties of the gas 
are determined by a weighted average over neighboring 
particles, formally discontinuous shocks are spread over 
a length determined by the smoothing length. Addition- 
ally, SPH relics on artificial viscosity to dissipate flows 
and produce the cor rect amount of e ntropy . Because 
of these restrictions, iPfrommer et al.l (|2006h developed 
a method that is able to identify shocks by examining 



the time-evolution of the entropy of individual SPH par- 
ticles. Comparing the instantaneous entropy injection 
rate to the characteristic time it takes a particle to cross 
the broadened shock surface, they are able to identify 
and calculate the instantaneous Mach number of shocks 
while remembering the pre-shock conditions. Therefore, 
the analysis can be performed on-the-fly and shock quan- 
tities were traced along with the usual hydrodynamic 
properties. These calculations use calibrations against 
"shock tube" simulations to derive the correct relation- 
ship between entropy injection rate and Mach number, 
which may vary with respect to different artificial viscos- 
ity implementations. 

To address all of the problems listed above, we have de- 
veloped a novel numerical algorithm capable of detecting 
and identifying shocks in the 3-D adaptive mesh refine- 
ment (AMR) grid-based code, Enzo. The use of AMR 
allows us to analyze unprecedented dynamic ranges with 
an advanced hydrodynamic code that is able to capture 
shocks exceedingly well. In this work we explore simu- 
lations with dynamic range of up to 2 16 = 65, 536, but 
we are not limited from going further in future work. 
Because of the complexity of the structure of AMR sim- 
ulations, it was necessary to develop several new numeri- 
cal algorithms to identify the shocks. This shock-finding 
analysis algorithm wil l be presented and compared to 
previ ous methods (e.g. IRvu et al.l [20031 IPfrommer et al.1 
2006). In order to validate and quantify the robustness 
of our method, we carry out a resolution study that in- 
cludes both mass and spatial resolutions that vary by 
factors of 16 and 64, respectively. 

We also propose a new method of characterizing shocks 
by their pre-shock overdensity and temperature. This 
then allows analysis that goes beyond the traditional in- 
ternal vs. external ( of filaments/cluste rs) shocks classi- 
fication suggested by IRvu et aTl (|2003l ) . By refining the 
temperature and density ranges examined, we are able to 
identify shocks that reside in voids, filaments, and halos. 
Additionally, by using temperature cuts, we can iden- 
tify the population of gas that is being shocked into the 
warm- hot intercluster medium (WHIM). 

After calculating the shock structure in a given sim- 
ulation, we are able to compute the amount of shock 
kinetic energy that is transferred to high-energy cosmic 
ray protons through diffusive shock acceleration. While 
the surface area of the large scale structure shocks are 
dominated by low pre-shock temperature and density, 
the bulk of the cosmic ray energy generation occurs in 
the centers of collapsed struc tures. Since stronger shock s 
will produce harder spectra (Blandfo rd fc Eichlerlll987| ). 
we expect that the strong accretion shocks could be the 
source of high energy cosmic rays. 

In Section 2 we describe the numerical methods used 
for both the cosmological simulations as well as the anal- 
ysis of the shock-finding proc ess. In Sec t ion 3 , we com- 
pare our algorithm to that of IRvu et alj (|2003l ) and test 
it using 3-D "shock tube" tests. Section 4 contains the 
main results of analyzing a large, (512 Mpc/h) 3 , cos- 
mological simulation with a peak spatial resolution of 
7.8 kpc/h. Section 5 describes the effects of spatial and 
mass resolution on the shock populations and cosmic ray 
acceleration. In Section 6, we discuss the limitations of 
our analysis, and in Section 7 we summarize our findings 
and discuss potential future directions. 
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2. METHODOLOGY 
2.1. The Enzo Code 

All simulations were run u si ng the Enzo cosmology 
code (iBrvan fc Norman1ll997allbl ;lN orman fc Brvanf l999: 
10 'Shea et al.l l2004f h While a full description can be 
found in the cited papers, we will review the key aspects 
that are of importance to this work. 

Enzo is a block-structured ad aptive mesh refinement 
(AMR; iBerger fc Colellal Il989f ) code that couples an 
Eulcrian hydrodynamics method that follows the gas 
dynamics with an N-Body particle mesh (PM) so lver 
(jEfstathiou et all Il985t iHocknev fc Eastwood! 1 19881) to 
follow the dark matter component. Enzo implements two 
hydrodynamic solvers. The first is a piecew ise parabolic 
method (PPM; IWoodward fc Colellal 1981 with cosmo- 
logical modifications bv lBrvan et al.l (|1995f l. The second 
is the method from the ZEUS magnetohydrodynamics 
code ( Stone fc Norman! 1992a. b). In this work we re- 
strict ourselves to the PPM method because of its supe- 
rior shock-capturing ability and lack of artificial viscosity. 

The AMR scheme within Enzo is handled by parti- 
tioning the simulation volume into 3D rectangular solid 
grids. Each of these grids contain a number of grid cells 
that set the spatial scale on which the hydrodynamics is 
solved. If a region of cells within a grid is determined to 
require higher resolution, as judged by a number of re- 
finement criteria including gas/dark matter overdensity, 
minimum resolution of the Jeans length, local gradients 
of density, pressure, or energy, shocks, or cooling time, 
then a minimum enclosing volume around those cells is 
created at the appropriate level of refinement. These 
newly created "child grids" can then themselves recur- 
sively become "parent grids" to yet another more highly 
refined region. This recursive nature does not set any 
restrictions to the number of grids or level of refinement. 
However, because structure formation leads to an enor- 
mous dynamic range we are limited by available compu- 
tational resources, a maximum level of refinement l max 
is instituted. 

In addition to being adaptive in space, Enzo imple- 
ments an adaptive time stepping algorithm. All grids on 
a given level are given advanced simultaneously with a 
maximum timestep such that the Courant condition is 
satisfied by all the cells on that level. This results in 
a hierarchy of timcsteps: a parent grid on level I is ad- 
vanced At(l), and then its subgrid(s) on level I + 1 are 
advance by one or more timesteps until they reach the 
same physical time as their parent grid. At this point, 
flux information is exchanged from child to parent grid 
in order to provide a more accurate solution to the hy- 
drodynamics on the parent grid. This procedure is done 
recursively until all grids are at the same physical time 
as the root grid, at which point the process is repeated 
until the stopping point of the calculation is reached. 

2.2. Shock-Finding Algorithm 

The bulk of our analysis relies on accurately identifying 
and quantifying the strength of shocks. After finding a 
shock, we would like to calculate its Mach number, which 
characterizes the strength of the shock. There are several 
methods that can be used in order to calculate the Mach 
number, including density, temper ature, velocity, or en- 
tropy jumps across the shock. As in lRvu et al.l(|2003l ). we 



use the Rankine-Hugoniot temperature jump conditions 
to calculate the Mach number. The temperature jump 
is preferable to density because it is more sensitive to 
Mach number, whereas the density jump quickly asymp- 
totes for strong shocks. The Mach number is solved for 

by 

T 2 (5A^ 2 -l)(A^ + 3) 

T x 16A4 2 ' 1 ' 

where T 2 and T\ are the post-shock (downstream) and 
pre-shock (upstream) temperatures, respectively. M. is 
the upstream Mach number. 

A cell is determined to have a shock if it meets the 
following requirements: 

V-w<0 (2) 
VT • VS > (3) 

T 2 > Tx (4) 

Pi > Pi: (5) 

where v is the velocity field, T is the temperature, p is the 
dens ity, and S = T/p 1 -1 is the entropy. In our analysis, 
as m IRvu et alJlpOOah . we have set a minimum prcshock 
temperature of T = 10 i K since the low-density gas in 
our cosmological simulations is assumed to be ionized (a 
reasonable assumption at z < 6). Therefore, any time 
the pre-shock temperature is lower than 10 4 K, the Mach 
number is calculated from the ratio of the post-shock 
temperature to 10 4 K. This introduction of a tempera- 
ture floor prevents us from drastically overestimating the 
accretion shock strength in adiabatic simulations. Future 
work will incorporate a self-consistent UV ionizing back- 
ground radiation. 

Now the task is to identify all of the shocks and 
their corresponding Mach numbers. The method that 
has been previo usly used in unigrid simulations (e.g. 
IRvu et ail l2003f ) is to loop through rows of cells along 
each of the coordinate axes and identify 1-D shock struc- 
tures in each direction. Contiguous cells that meet the 
requirements above in Eqs. 2-5 are then combined into 
a single shock structure with the cell of maximum con- 
vergence marked at the center. The pre- and post-shock 
cells are identified as first cells outside of the shock struc- 
ture. If a center is marked in more than one direction, 
the maximum calculated Mach number of the three pos- 
sible is taken to be the true Mach number. Because of 
this, we would expect errors to arise when examining 
shocks whose direction of motion is not oriented along a 
coordinate axis. To address this issue, we have designed 
an algorithm that does not suffer from this limitation. 

In our method, we first determine the direction of 
shock propagation from the local temperature gradients, 
making the assumption that the shock-induced temper- 
ature gradient overwhelms the background temperature 
gradient. We then search the cells along the tempera- 
ture gradient for the pre- and post-shock cells. If we 
find a neighboring cell to have a more convergent flow, 
that cell is marked as the center and we move outwards 
from it. This guarantees that the analysis is anchored 
to the center of the shock. Once the furthest pre- and 
post-shock cells are found, the temperatures are taken 
and the Mach number is calculated from Equation 1. A 
two-dimensional analog of this process is shown in Fig- 
ure [T] Because our algorithm is not confined to operate 
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Fig. 1. — A 2-D cartoon of our shock analysis algorithm. The 
shock centers are shown as dark blue cells, while the pre- and post- 
shock cells are outlined in thick black. The AMR resolution level 
is seen by varying grid-cell sizes. 

along coordinate axes, the calculated Mach numbers pro- 
vide a more accurat e description of th e underlying shock 
properties than the lRvu et a l. ( 2003) method. 

Specifically, in situations where there are weak shocks 
or complex flow velocities, using the coordinate-split 
approach may allow for an excess of shocks since the 
direction of the Shockwave is not taken into account. 
Our method picks a single direction that a shock could 
be propagating, given a specific temperature gradient. 
I Rvu et all (|2003D claim that their shocks are spread out 
over 2-3 cells, of which one is marked as the center. For 
the other 1-2 cells, the coordinate split approach may 
mis-identify these other 1-2 cells as low Mach number 
shocks due to normal temperature gradients. This would 
lead to an over-prediction of low Mach number shocks. 

This process is further complicated by the use of adap- 
tive mesh refinement (AMR), primarily because neigh- 
boring cells are not necessarily at the same level of re- 
finement. This occurs most often at the site of accretion 
shocks onto halos and filaments where the density gra- 
dient, upon which the refinement criteria are based, is 
largest. Therefore, knowledge of the grid hierarchy must 
be used. We incorporate this into our algorithm and al- 
low for a neighboring cell to be any cell at the same or 
lower level of refinement. We do not allow the algorithm 
to search for neighbors at higher levels (smaller grid cells) 
since one cell will have multiple neighbors. If this situ- 
ation occurs, we use the neighboring cell on the same 
level. Because of this requirement, we must perform our 
analysis on the most highly refined grids first, and move 
to progressively coarser levels of resolution. 

2.3. Cosmic Ray Acceleration Models 

Following the method proposed in IRvu et alj (|2003) , 
we now seek to determine the amount of kinetic energy 
that is converted into heating of the gas and accelerating 
cosmic rays. We begin with calculating the total kinetic 
energy flux through a shock surface. The kinetic energy 
flux associated with a mass flux of p\M.c s is: 

1 
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Fig. 2. — Fractional efficiency of gas thermalizat i on an d cosmic 
ray acceleration, from the models by Kane: & Jones (2007). 5o(A4) 
is the gas thermalization fraction expected from the Rankine- 
Hugoniot jump conditions. 8(J\4) and r){M) are the gas thermal- 
ization fraction and cosmic ray accelerati on frac t ion, re spectively, 
from the non-linear calculations of Kang & Jones (2007). The two 
panels show results assuming different compositions of the pre- 
shock plasma. Top: thermal gas with no cosmic ray population. 
Bottom: thermal gas with a pre-exisiting cosmic ray population 
having PqrI Pgas = 0.3, where Pqr and P g are the cosmic ray 
and total gas pressure, respectively. 

where p\ is the pre-shock density and c s is the sounds 
speed in the pre-shock gas. From this total incoming ki- 
netic energy flux, a fraction will be used in the thermal- 
ization of the g as and the acceler ation of cosmic rays. 
In keeping with I Rvu et all (|2003h . we will denote the 
amount of energy per unit time used to heat the gas 
and accelerate cosmic rays as Jte and fcR, respectively. 
In the case of a purely hydrodynamical shock without 
the inclusion of cosmic ray feedback, the fractional ther- 
malization Sq(M) can be determined by the Rankinc- 
Hugoniot jump conditions, 



S Q (M) = 



e-TE.2 — CTE.l 



''2 



(7) 



KE 



- Pl {Mc,y 



(6) 



where &te,i and cte.i are the thermal energy densities 
in the pre and post-shock regions, respectively. 

With the inclusion of cosmic rays, there is no simple 
analytical form for the fractional thermalization of the 
gas, which depends on magnetic field orientation, tur- 
bulence, and the pre-shock cosmic ray population. In- 
stead, we adopt the results of ID diffusive shock ac - 
celeration (DSA) simulations by I Kang fc J ones! (|2007l) . 
The time-asymptotic values of the fractional thermal- 
ization, 5(M) = Jte/Jke, and fractional CR acceler- 
ation, ri{M) = JcrI Ike, were found to be nearly self- 
similar for the temperatures and shock velocities consid- 
ered. These simulations also accounted for whether or 
not the pre-shock medium had preexisting cosmic rays. 
With a preexisting cosmic ray population, the fractional 
energy deposited into CRs increases dramatically at low 
Mach numbers because it is much easier to accelerate an 
existing power-law distribution of particles than a ther- 
mal distribution of particles. Shown in Figure [21 are the 
results of the Kang & Jones DSA simulations for a pop- 
ulation with no preexisting cosmic rays and one in which 



5 



TABLE 1 
Simulation Parameters 



Name 


Lbox 




Imax 








Q m 


erg 


rj/ul024 


100 


97.7 





5.877 x 10 7 


97.7 


0.043 


0.27 


0.8 


SF Light Cone 


512 


1 


7 


7.228 x 10 10 


7.8 


0.04 


0.3 


0.9 


mlJ8 


256 


1 


8 


6.224 x 10 10 


3.9 


0.0441 


0.268 


0.9 


mlJ6 


256 


1 


6 


6.224 x 10 1() 


15.6 


0.0441 


0.268 


0.9 


rolJ4 


256 


1 


i 


6.224 x 10 10 


62.4 


0.0441 


0.268 


0.9 


ro4J8 


256 


500 


8 


7.781 x 10 9 


3.9 


0.0441 


0.268 


0.9 


m,4J6 


256 


500 


6 


7.781 x 10 9 


15.6 


0.0441 


0.268 


0.9 


to4J4 


256 


500 


4 


7.781 x 10 9 


62.4 


0.0441 


0.268 


0.9 


ml6J8 


256 


250 


8 


9.726 x 10 8 


3.9 


0.0441 


0.268 


0.9 


ml6J6 


256 


250 


6 


9.726 x 10 8 


15.6 


0.0441 


0.268 


0.9 


ml6J4 


256 


250 


i 


9.726 x 10 8 


62.4 


0.0441 


0.268 


0.9 



Note. — Lb ox is the simulation box size in comoving Mpc/h. J\xrq is the effective root grid resolution (the m4 and ml6 series of 
calculations use one and two static nested grids, respectively). l max is the maximum level of AMR allowed in the simulation. M dm is the 
dark matter particle mass (in the static nested grids for the m4 and ml6 series of runs) in M^/h. Ax ma x is the peak spatial resolution 
in comoving kpc/h. Cli, and f2 m are the fractional densities of baryons and matter compared to the critical density (Q\ = 1 — Qm in all 
simulations, so Qq = 1). as is the power spectrum normalization of the mass fluctuation in a comoving 8 Mpc sphere. 



CRs existed initially with Pcn/Pg ~ 0.3, where Pen and 
P g are the cosmic ray and total gas pressure, respectively. 
The sum of 6(A4), ij(Ai), and the remaining fraction of 
kinetic energy in the gas is equal to one, conserving en- 
ergy. 

We do not track the cosmic ray population in our sim- 
ulations at present, and as a result we are unable to con- 
strain the amount of preexisting CRs in the pre-shock 
medium. Therefore, we can think of our results from the 
two scenarios shown in Figure [2] as bracketing the likely 
range of efficiencies. Additionally, these efficiency mod- 
els are only valid for situations where the shock normal is 
parallel to the magnetic field. Any deviation from these 
ideal conditions will likely reduce the e fficiency of cosmic 
ray acceleration (|Kang fc Jo nes 20071 ) . so one can view 
the results described later in this paper as upper limits 
on cosmic ray injection efficiency. 

2.4. Simulations 

We constructed three distinct sets of cosmological sim- 
ulations for this project. A summary of some of the 
simulation parameters is given in Table [T] First, we have 
a simulation that was devised as an analog of the uni- 
grid numerical simulation by Ryu et al. (2003). For this 
simulation , we us ed identical cosmological parameters to 
iRvu et alJ (|2003f) in order to provide a reference simula- 
tion to compare our new shock-finding method with pre- 
vious work. The cosmological parameters for this sim- 
ulation, ryul024, are: Q, B m = 0.043, Q DM = 0.227, 
n A = 0.73, h = ff /(100 km s^Mpc' 1 ) = 0.7, and 
er 8 = 0.8, which are broadly con sistent with WMAP Year 
5 results (jKomatsu et al.l)2008[ ). The comoving size of 
the simulation volume is (100 Mpc/h) 3 and is discretized 
into 1024 3 cells, giving a comoving spatial resolution of 
97. 7kpc/h. It also employs 512 3 dark matter particles 
with a 1024 3 grid. In order to reproduce the Ryu et al. re- 
sults as closely possible, we did not use AMR techniques. 
The simulation was initialized with an Eisenstein & Hu 
(1999) power spectrum with a spectral index of n = 1.0 
at z=99 and the simulation states were output in 20 times 
between z = 10 and z — 0. The analysis of this simula- 
tion is described in Section 3.2. 

Our main results in this work focus on the analysis of 
a (512 Mpc/h) 3 volume that utilizes a 512 3 root grid 



and up to 7 levels of AMR. It is referred to as the 
"Santa Fe Light Cone," a nd ha s been previously de- 
scribed by lHallman et al.l (2007) . This simulation has 
a peak spatial resolution of 7.8kpc/h and a dynamic 
range of 65, 536. The cosmological parameters used were: 

n M = 0.3, n BM = 0.04, n C DM = o.26 , n A = 0.7, 

h = iZ /(100 km s^Mpc' 1 ) = 0.7, and a 8 = 0.9, and 
employs a Eisenstein & Hu (1999) power spectrum with 
a spectral index n — 1.0. Cells are refined whenever the 
baryon or dark matter density increased by a factor of 8 
beyond the previous level. Because the simulation then 
refines by a factor of 8 in volume, the average mass per 
grid cell stays roughly constant. The simulation was ini- 
tialized at z = 99 and was run to z = 0. The analysis of 
this simulation is described in Section 4. 

In order to study the effects of spatial and dark 
matter mass resolution, we have performed a suite 
of simulations that vary these factors, and illustrate 
their results in Section 5. These simulations are known 
as "nested grid" simulations. An initial cosmological 
simulation is run at low resolution. The most massive 
halo at z = is found and the simulation is re-centered 
at the final location. The simulation is then re-run, 
while only adaptively refining a region that bounds all 
dark matter particles that eventually are inside the most 
massive halo. Therefore, the focus of the simulation is 
only on the inner portion of the initial volume. With 
this capability, we are able to modify the root grid 
and peak spatial resolution for this subvolume and 
study their direct effects on the evolution of a single 
cluster. In our simulations, we initialize a (256 Mpc/h) 3 
volume with 256 3 root grid cells. From that, we only 
adaptively refine in a (32 Mpc/h) 3 subvolume. Within 
the subvolume, we add up to two static nested grids, 
with more highly refined dark matter particles and 
gas cells. A list of all simulation parameters used is 
given in Table [TJ The cosmological parameters used 
are: Q M = 0.268, Q BM = 0.0441, n CD M = 0.2239 , 
ft A = 0.732, h = H /(100 km s^Mpc" 1 ) = 0.704, and 
<t$ = 0.9. These parameters are the WMAP year 3 
parameters (jSpergel et alJ 12003) but with a somewhat 
higher erg. 
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3. VALIDATION OF SHOCK-FINDING METHOD 
3.1. Shock Tube Test 

In order to verify that our shock-finding algorithm is 
accurate, we have performed a suite of 3-D AMR shock 
tube tests. In these tests we have varied the Mach 
number as well as orientation with respect to the co- 
ordinate axis. The setup of this test pro blem is de- 
scribed in lMihalas fc Weibel Mihalasl (|1984l) . It consists 
of a stationary, uniform pre-shock medium. The shock 
is then introduced via boundary conditions that match 
the Rankine-Hugoniot Jump conditions for a given Mach 
number. The volume is then allowed to adaptively refine 
up to 2 levels, using shocks as a criteria for refinement. 
We have chosen to adaptively refine based on shock lo- 
cations (i.e. strong pressure jumps) instead of density 
because this should introduce a complicated AMR topol- 
ogy in order to test the robustness of our algorithm. This 
forces us to traverse different levels of refinement for pre- 
and post-shock quantities. 

In order to change the direction of shock propaga- 
tion, we change the time at which a given boundary cell 
changes from uniform to "shocked." Using this proce- 
dure, we vary the shock propagation vector over both 9 
and <j), which are angles off of the x-z and x-y planes, 
respectively. In addition to the three on-axis scenar- 
ios, we vary 9 and </> over all permutations of the angles 
0, 7r/8, 7r/6, and7r/4. For each shock propagation direc- 
tion, we then vary the input Mach number over M. =2, 
5, 30, and 100. 

The general result from this study is that our shock- 
finding algorithm is very accurate. As shown in Figure 
[31 if we make a histogram of the ratio of calculated Mach 
number to expected Mach number, and normalize it so 
that the area under the curve is equal to 1, the result is 
both accurate and precise. In Figure [31 we created the 
histogram by summing over all orientations of the shock 
of a given Mach number. As one can see, the peak is 
centered around 1.0, with an average sample standard 
deviation of less than 0.06. We have examined the aver- 
age Mach number and standard deviation as a function 
of angle and have found no discernible trend or bias. 

Additionally, due to the manner in which we set up the 
propagating shock, small inhomogeneities arise that are 
likely the cause of much of the calculated scatter. This 
is because we introduce the shock from the boundary 
conditions which do not explicitly keep the leading edge 
of the shock as a perfect discontinuity. Therefore the 
accuracy of our shock finding algorithm is likely better 
than that shown in Figure [3] 

3.2. Comparison to \Rvu et al\ II 200 A ) 

I n order to t e st our analysis against previous work done 
by IRvu et alj (|2003ft . we generated a 1024 3 fixed grid 
simulation with identical cosmological parameters and 
spatial resolution as their most highly-resolved calcula- 
tion. This simulation is dexcribed in detail in Section 
2.4. We expect to see a difference in results from the 
shock-finding method and from underlying differences in 
the hydrodynamical solvers. Enzo uses the Piecewise 
Parabolic Met hod, which captur es shocks across a single 
zone, whereas IRvu et all ((2003) use the total variation 
diminishing (TVD) method, which spreads shocks over 
approximately two cells. 
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Fig. 3. — Distributions of the ratio of calculated Mach number 
to expected mach number from off-axis 3D AMR shock tube test 
problems. Shock surface area, S, distributions were averaged over 
all orientations for the each individual Mach number, and normal- 
ized so that the area under the curve is 1. Varying lines correspond 
to Mach numbers of 2 (dash-dotted), 5 (dashed), 30 (dotted), and 
100 (solid). Sample standard deviations from 1.0 are all less than 
0.06 

In order to study the diff erences between our shock- 
finding methods and those of lRvu et alj (120031) ■ w e mimic 
the top half of Figure 5 from IRvu et all (|2003h m our 
Figured! However, in addition to using our new shock- 
finding algorithm that searches along temperature gradi- 
ents, we in clude a coo r dinate -split analysis that is similar 
to that of IRvu et all ((2003). As we claimed in Section 
2.2, using a coordinate-split approach overpredicts the 
number of low Mach number and internal shocks. For ex- 
ternal, low Mach number shocks, the difference is roughly 
a factor of 3, which agrees with our hypothesis that the 
coordinate-split approach identifies cells that are associ- 
ated with a strong shock to the center of a weak shock. 
The difference in the internal shocks spans all Mach num- 
bers because the flow to be very complex, making it easy 
to mistake a normal temperature gradient with that of 
a shock. We have also studied the integrated kinetic 
flux through sh ock surface s and find it to be in general 
agreement with IRvu et all (|2003l ) . We will study this in 
more detail in the future when we include this analysis 
"on-the-fly." 

While the shock surface area distributions are good 
indicators of qualitative differences, we now quantify 
these resu l ts. T his is done by recreating Table 1 from 
IRvu et all (|2003h in our Tabled For this portion of the 
analysis, we mimic the Mach number floor requirement 
that M > 1.5, whi ch was used to re duce the effects of 
complex flow in the IRvu et all (|2003h analysis. First, it 
is instructive to give a physical motivation for these pa- 
rameters. The quantity 1/5* can be thought of as a mean 
separation of shocks because it is the simulation volume 
divided by the total shock surface area. This gives it 
units of comoving Mpc. The ratio of external and inter- 
nal shocks gives the reader an intuition as to where the 
majority of the shocks are occuring. Note that as the 
redshift decreases, the relative amount of internal shocks 
increases, indicative of the increase in shocks within ha- 
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shocks for z = in the ryul024 simulation. The two methods of shock finding, coordinate split (dotted line) and coordinate unsplit (solid 
line), are shown. At low Mach numbers for external shocks and for all i ntern al shocks we see a significant overprediction in the number of 
shocks when using the coordinate split method described in lRvu et a l. (2003) 



TABLE 2 
Mean Shock Quantities 



z 


1/S 




(M ext ) 


(M int ) 


1/ S ex t 


1/ Si n t 


0.0 


6.235 


3.550 


12.65 


3.767 


7.992 


28.37 


0.25 


6.519 


4.312 


13.02 


3.961 


8.030 


34.63 


0.50 


6.886 


5.196 


12.86 


4.119 


8.211 


42.67 


0.75 


7.301 


6.248 


12.61 


4.225 


8.470 


52.92 


1.0 


7.767 


7.442 


12.26 


4.310 


8.811 


65.57 


1.25 


8.297 


8.743 


11.81 


4.399 


9.246 


80.84 


1.50 


8.884 


10.18 


11.34 


4.412 


9.756 


99.38 


1.75 


9.546 


11.59 


10.89 


4.411 


10.37 


120.2 


2.0 


10.31 


13.04 


10.50 


7.679 


11.10 


144.8 



Note. — Mean shock quantities, z is the redshift of the simu- 
lation. 1/S is the mean comoving length between shock surfaces in 
units of h~ 1 Mpc. S e xt/Si n t is the ratio of shock surface area for 
external to internal shocks. {M ex t) and (Mint) are the surface area- 
weighted mean of the external and internal shock Mach number, 
respectively. 1/Sext and 1/Si n t are the average comoving distance 
between external and internal shocks, respectively, in h~ 1 Mpc. 



los and the measured amount of matter in large halos. 
The average quantities are surface area weighted means 
of the quantity in question. A subscript of ext or int 
denotes that only external or internal shocks were used, 
respectively. External shocks are those with pre-shock 
temperatures less than 10 4 K while internal shocks are 
those with pre-shock temperature greater than 10 4 K . 

In comparing our Table [2 to Table 1 in Ryu 
et al. (2003), we find that we predict a higher average 
Mach number and higher mean comoving distances be- 
tween shocks for both internal and external shocks. The 
ratio between our average external Mach number and 
that found in Ryu et al. ranges between 1.54 and 1.6, 
while that of the internal Mach number (disregarding 
z = 2.0) ranges between 1.3 and 1.5. We have disre- 
garded z = 2.0 because there is a large amount of merg- 
ing between z — 2.0 and 1.75, significantly raising the 
internal temperature of many of the large clusters, in- 
creasing the sound speed and decreasing the Mach num- 
ber. Therefore, we believe that our particular realization 



of this volume had later mergers than that of iRyu et al.l 
(l2003h . 

The differences in the average Mach numbers as well as 
the increase in mean comoving distance between shocks 
is almost entirely due to the use of a coordinate split al- 
gorithm vs. a coordinate unsplit algorithm. The identifi- 
cation of many more low Mach number shocks increases 
the frequency, thus decreasing the comoving length be- 
tween shocks. Therefore for future studies, this difference 
must be taken into account. 

4. RESULTS FOR THE SANTA FE LIGHT CONE VOLUME 

Now that we have outlined our improved shock finding 
algorithm, we apply it to a large cosmological simulation 
encompassing a volume of (512 Mpc/h) 3 . This simula- 
tion, called the "Santa Fe Light Cone," was described 
previously by Hallman et al. (2007). This represents the 
first time that a large cosmological volume with superb 
spatial resolution has been studied for its shock and cos- 
mic ray properties. Whereas previous studies were only 
able to study a small nu mber (~ 10) of cluster s due to a 
small cosmological box (Pfromm er et al.ll2006| ). we have 
over 9000 halos with M ha i > 5 x 10 13 M Q , and over 200 
with Mhaio > 5 x 10 14 M Q . This allows us to perform a 
statistical study of cosmological shocks unlike any that 
has been done previously. Both the increase in volume 
(by a factor of ~ 125) and an enhanced spatial reso- 
lution over previous unigrid/SPH simulations allow un- 
precedented detail in our calculations. 

We begin by outlining the shock distribution and how 
it can be thought of as a new way to view large scale 
structure formation in the Universe. We do this by break- 
ing the distributions down by temperature and density 
cuts, which further illuminates the underlying dynamics. 
From there, we apply the DSA cosmic ray acceleration 
model and determine what phase of gas will contribute 
most to the acceleration of cosmic rays. Finally, we es- 
timate the global fraction of kinetic energy that is pro- 
cessed through shocks that is devoted to the acceleration 
of cosmic rays in an effort to determine their possible 
dynamical effects on gas behavior in galaxy clusters. 



Fig. 5. — Projections of a 2.8 X 10 h Mq cluster from the "Santa Fe Lightcone." Mach number (top-left) is weighted by the injected 
cosmic ray flux. Injected cosmic ray flux (top-right) is in units of ergs/(s h~ 2 Mpc~ 2 ). Baryon column density (bottom-left) is in units of 
Mq / (h~ 2 Mpc 2 ) . Mass-weighted temperature (bottom-right) is in units of Kelvin. The total size of the projected volume is (32 /i -1 Mpc) 3 . 
All panels show logarithmic quantities. 



4.1. Shock Frequencies 

As was done in lRvu et al.1 (|2003f ). we calculate the sur- 
face area of all shocks in a given logarithmic Mach num- 
ber interval. However, instead of only classifying shocks 
as internal or external depending on their prcshock tem- 
perature, we break the distribution into logarithmic tem- 
perature and density cuts that can be postprocessed to 
examine any subset of the p or T phase space for the 
entire computational volume. Primarily, we create sev- 
eral physically motivated temperature and density cuts, 
which are outlined in Table [3] Note that the gas in the 
T < lO^K heading is artificial since we do not include a 
UV background. This temperature range traces gas that 
has not been previously shock heated. In addition to 
studying the physical properties of the preshock region, 
we study the evolution of the distributions as a function 
of redshift . 

Figure [5] shows a projection of the Mach number for the 
largest cluster in the simulation (2.8x 10 15 M Q ), weighted 
by cosmic ray acceleration rate. This allows us to see 
the structure of cosmological shocks. By weighting the 



TABLE 3 
Temperature-Density Phase Space 



Location 


Temperature Range 


Overdensity Range 


Voids 


T < 10 4 K 


S< 1 


Filaments 


10 4 if < T < 10 6 K 


1 < S < 100 


Clusters 


10 6 K < T < 10 S K 


100 < 8 < 10 3 


Cluster Cores 


T > 10 a K 


S > 10 3 



Note. — Approximate ranges for pre-shock temperature 
T or pre-shock overdensity <5 = ph/{pb) f° r general large scale 
structures. 



projection by cosmic ray acceleration rate, we see both 
the external high-Mach number shocks and the internal 
shocks, since the internal shocks' weights are higher. In 
the other three panels, we show the injected cosmic ray 
flux, density, and mass-weighted temperature. 

4.1.1. Density & Temperature Ranges 

We now ex pand the c l assifi cation of external and in- 
ternal shocks iRvu et all (|2003f ) by examining the shock 
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Fig. 6. — DifTcrential shock surface area as a function of logarithmic mach number bins for varying pre-shock gas phases. Pre-shock gas 
overdensity (top) is divided into several ranges that differentiate the overall distribution. Pre-shock gas temperature (bottom) differentiates 
the different types of shocks (i.e. accretion, merger). Both distributions are shown for three redshifts: z = 3 (left), z = 1.5 (middle) and 
z = (right). 



Mach number distributions in varying temperature and 
density ranges. This will provide a more complete de- 
scription of where these shocks arise in structure for- 
mation than in previously published analyses. In Fig- 
ure [51 the shock frequency is plotted for a range of 
density and temperature cuts. At z = 3, we see that 
shock surface area distribution is dominated by shocks 
with low temperature/low overdensity pre-shock quanti- 
ties. These represent the accretion shocks onto filaments 
and proto-clusters. As the simulation evolves, the dis- 
tribution becomes bimodal with components from both 
low pre-shock temperature, high-Mach number accretion 
shocks and high-temperature, low-Mach number merger 
shocks. 

The temperature cuts each have a characteristic Mach 
number cutoff that increases with decreasing tempera- 
ture. This cutoff is due to the maximum temperature 
jump that is possible with a given pre-shock temper- 
ature. Therefore, since the maximum temperature in 
the simulation is ~ 10 8 iv~ (determined by the mass of the 
largest cluster), a temperature jump from 10 e K will re- 
sult in a A4 « 18, very close to the cutoff seen at z = 
for 10 6 K < T < 10 7 K. Similarly, the cutoffs for lower 
pre-shock temperatures indicate the largest temperature 
jumps for each population. At higher redshifts, these 
temperature cutoffs decrease due to the lower maximum 
temperatures present in the simulation. Therefore, the 
movement of these cutoffs tell us about the temperature 
evolution of the simulation. 

Additionally, the Mach number associated with the 



peaks in the shock frequency distribution can be used 
to determine the mathematical mode of the post-shock 
temperature distribution. For z = 0, these peaks cor- 
respond to post-shock regions with ~fewxl0 6 -ftT for 
pre-shock temperatures T\ < 10 6 K . Therefore the ma- 
jority of these shocks are heating the pre-shock gas to 
WHIM temperatures in filaments. For Ti > 10 6 A', the 
peak Mach numbers correspond to post-shock tempera- 
tures of T 2 « 2 x 10 7 — 1.5 x 10 8 K. These are complex 
flow and subhalo merger shocks in the interior regions of 
clusters. 

If we instead examine the varying density cuts, similar 
results are observed. At high redshifts, we see that the 
dominating accretion shocks (high Mach number shocks) 
have pre-shock overdensities of S ~ 1 — 10. This is be- 
cause of the relative paucity of large-scale halos and fil- 
aments and, thus, relatively shallow gravitational poten- 
tial wells. The infalling gas will get much closer to the 
accretor and therefore denser before shocking. As we 
move to lower redshift, the S < 1 shocks begin to domi- 
nate because we are shocking further out into the voids. 

For the interior cluster shocks, there are three regimes 
that present themselves in the analysis. If we examine 
z = 3 with 10 < S < 100, there are plateaus near M. « 
2-4 and M « 10 - 70. It is difficult to determine what 
the post-shock density will be because of the insensitivity 
of the density contrast at high Mach numbers (p2 / pi — » 4 
for M. » 3). However, it is likely that the two high 
Mach number shock plateaus correspond to filaments for 
M. « 10 and clusters at the virial radius for M ps 70. 
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Fig. 7. — Comoving shock surface area normalized by the sim- 
ulation volume as a function of Mach number with varying red- 
shifts. Three regions are suggested corresponding to internal clus- 
ter merger shocks (red), accretion shocks onto filaments (green), 
and accretion shocks onto clusters (blue). 

The low Mach number shocks are most likely interior flow 
shocks. 

At late times, all of the intermediate pre-shock density 
regions have bimodal distributions. The high Mach peak 
corresponds to density contrasts of 4, while the low M 
corresponds to jumps of ~ 2. Therefore, we are likely 
looking at merger and complex flow shocks, respectively. 

4.1.2. Redshift Evolution of Shock Properties 

There are three primary populations of shocks that 
we see evolve through time, as seen in Figure [7] There 
are accretion shocks onto clusters, accretion shocks onto 
filaments, and merger and complex flow shocks within 
clusters and filaments. These are outlined by the blue, 
green, and red shadings in Figure [JJ and their qualita- 
tive behavior can give useful insights as to the evolution 
of large scale structure. Let us analyze each of these 
populations separately. To determine the origin of these 
populations, we have examined slices and projections of 
the data and compared the Mach number of the cell to 
it's location with respect to large scale structure. 

First, at early times we see a small peak at very high 
Mach numbers that denotes shocks onto collapsing ha- 
los. This corresponds to gas that has previously been un- 
touched by shocks falling directly onto the proto-cluster 
gas, with temperature jumps from hundreds of Kelvin 
to lO 6 ^ (Note the Mach numbers are still calculated 
with a temperature floor of 10 4 K). We see that as the 
universe evolves, the strongest shocks in the simulation 
become stronger. This is due to the mass of the clus- 
ters increasing with time, providing a larger gravitational 
force pulling the material onto the halo. We also see that 
this peak increases in shock frequency while slowly mov- 
ing to slightly lower Mach numbers. Because the mass 
function cuts off exponentially at high mass, the number 
of small halos heavily outweighs the large halos. These 
smaller halos have lower free-fall speeds at the radius of 
the accretion shock, leading to a smaller Mach number. 
Therefore the large number of weaker shocks dominate 
the net surface area distribution. 
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Fig. 8. — Redshift evolution of the amount of kinetic energy 
processed by shocks as a function of Mach number. Redshift de- 
creases from z = 3 (red) to z = (black). The decrease in flux 
at late times for M > 10 signals the epoch at which dark energy 
becomes dominant. 

Second, the shocks onto filaments begin at Mach num- 
bers of M. ~ 6 and move towards M ~ 20 at late times. 
The surface area of these shocks are much larger at early 
times because the surface area of a cylinder per unit vol- 
ume is larger than that of a sphere as well as an increased 
number of filaments with respect to halos (there are sev- 
eral filaments that feed into a single halo). The strength 
of these shocks grow with the increase in size of the fila- 
ments. 

Finally, the low Mach number shocks (A4 < 3) due to 
halo mergers and complex flow are nearly non-existent at 
high redshifts. However, as large halos collapse and start 
to merge, the shock surface area also increases. There- 
fore this evolution traces the strength and frequency of 
merger shocks. 

4.2. Cosmic Ray Energy Injection 

The thermal history of the large scale structure in the 
Universe is primarily determined by the conversion of 
gravitational potential energy into kinetic energy, which 
is subsequently converted to heating gas and the acceler- 
ation of cosmic rays. Here we present results of our ap- 
plication of the cosmic ray acceleration model described 
in Section 2.3 to the "Santa Fe Light Cone." 

4.2.1. Function of Redshift 

T he first result is that, as in prev i ous studies 
(e.glRvu et alJ 120031 : IPfrommer et "all 120061 : iKang et all 
l2007fl . the most important Mach number shocks in terms 
of cosmic ray acceleration are at A4 w 2 — 4. This may 
seem surprising given that the surface area of shocks is 
dominated by high Mach shocks. However, the amount of 
energy dissipation is the product of the mass flux through 
the shocks and the Mach number. The large accretion 
shocks at early times (z < 3) have already consumed a 
large fraction of the gas in voids. This leaves very lit- 
tle mass at low densities to be processed by the most 
massive halos. This is in contrast to the low Mach com- 
plex flow shocks within the clusters. These process very 
large amounts of mass and kinetic energy, and therefore 
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Fig. 9. — Out of the incoming total kinetic energy of the shocks 
(solid line), the relative amount of energy devoted to the acceler- 
ation of cosmic rays for both models with (dashed line) and with- 
out (dotted line) a pr e-existing CR population, as predicted by the 
Kang & Jones (2007) diffusive shock acceleration model. 

experience very high thermalization and acceleration of 
cosmic rays even with lower efficiency. Cosmic rays from 
these low Mach number shocks will, however, have a 
steep energy spectrum and dissipate their energy rela- 
tively quickly comp ared to strong accretion shocks (e.g. 
iMimati et aLll200l( ). 

Figure [5] shows a distribution function of the kinetic en- 
ergy processed through shocks per comoving (Mpc/h) 3 
as a function of redshift where the height of the distri- 
bution function gives the differential amount of kinetic 
energy processed by shock for a given Mach number bin. 
As the simulation evolves to z = 0.5, there is a mono- 
tonic increase in the average kinetic energy density pro- 
cessed. Both the low-Mach complex flow and high-Mach 
accretion shocks increase by factors of 10 — 100. This 
monotonic increase stops at z ~ 0.5 because of the dom- 
inance of dark energy in a AC DM universe at this epoch, 
resulting in a decreased merger of, and accretion onto, 
the highest-mass halos. Therefore, the number of ac- 
cretion shocks characterized by high Mach numbers will 
decrease. Compounding this effect is the slow evacuation 
of the voids and the lack of additional mass to accrete. 

By applying the diffusive shock acceleration model, we 
can estimate how much of this energy is put into gas ther- 
malization versus the acceleration of cosmic rays within 
the confines of the model. This acts as a first estimate of 
the energy injection into cosmic rays, and should not be 
taken as the final on the subject. Cosmic ray injection is 
a highly non-linear process that is not fully understood. 
Further work on this model is needed. 

Figure [HI shows the relative amounts of energy dissi- 
pated for the two different models involving either no 
pre-existing cosmic rays or an initial amount of cosmic 
rays such that Pcn/Pg ~ 0.3. As one can see, the rela- 
tive amount of cosmic ray acceleration vs. thermalization 
heavily depends on the assumed inputs of the underlying 
DSA model. Until we are able to track the cosmic ray 
pressure within our simulations, we are resigned to give 
these rough limits of cosmic ray acceleration. 



4.2.2. Variation of Cosmic Ray Injection Efficiencies With 
Gas Properties 

Separation of distribution functions showing thermal- 
ization as a function of both temperature and density 
provides valuable insight into the physical processes oc- 
curring in the simulation. In Figure flOl we see that there 
are two primary modes of kinetic energy flux at z = 3. 
For M. < 2, the thermalization is dominated by shocks 
at 100 < 5 < 10 4 and T > 10 6 K. These shocks are likely 
within the largest filaments and the first clusters. At 
higher Mach numbers, A4 > 6 — 7, the thermalization is 
dominated by gas at T < 10 6 K and 8 ~ 10 - 100. This 
points towards accretion shocks onto filaments and the 
heating of the WHIM. If we use the peaks in each temper- 
ature cut up to T ~ 10 6 if to estimate the Mach number, 
we can calculate the post-shock temperature for these 
shocks to be 1 — 3 x 10 6 K. This reinforces the thought 
that these shocks are heating the WHIM. Shocks in this 
range of Mach numbers are also seen in Figure [5] as sur- 
rounding the filaments. 

At later times, the entire distribution shifts to higher 
thermalization rates due to the collapse of structures. 
Low Mach numbers are again dominated by complex 
flows within clusters. By examining the shocks with pre- 
shock temperature of less than 1Q 5 K as well as the red- 
shift evolution from Figure EJ we are able to verify that 
the high Mach number accretion shocks are becoming 
less important due to the separation of the voids from 
the clusters after z ~ 0.5. In the overdensity cut that 
corresponds to 1 < S < 10, we see a shift from a peak 
at high Mach numbers to small Mach numbers as the 
relative importance of accretion and me rgers switch. 

If w e compare our results to those of iPfrommer et aTl 
we see a good agreement at low Mach numbers. 
IPfrommer et al.l (|2006l ) found shocks as strong as M. ~ 
10 3 . However, we never see shocks above M. « 200. 
This is likely due to the lack of a temperature floor in 
their simulation, which thus allows a higher numerical 
value for the Mach number. These shocks are likely not 
realized in the real universe due to the presence of a 
ubiquitous ionizing radiation background that will keep 
gas above 10 4 K. 

Finally, we can compare our results to recent work by 
Kang et al. (2 007), who u sed a unigrid calculation sim- 
ilar to that of iRvu et a l. (2003), but included radiative 
processes, star formation, and a relaxed minimum tem- 
perature floor. Again, this relaxation of the temperature 
floor to (in their case) the CMB temperature resulted in 
very high Mach numbers - up to M. > 10 4 . This corre- 
sponds to a temperature jump by a factor of ~ 3 x 10 7 , 
a jump from 3K to 10 S K (the minimum and maximum 
temperatures in the simulation). At low Mach numbers, 
our re sults are very similar to those of iKang fc Jonesl 

(pool . 

5. EFFECT OF MASS AND SPATIAL RESOLUTION ON 
COSMIC RAY ACCELERATION EFFICIENCY 

In order to quantify the robustness of our simulations 
with respect to mass and spatial resolution, we perform 
a series of simulations where the mass and spatial reso- 
lution of a single galaxy cluster are varied over a wide 
range of parameter space. Two parameters are varied 
in this study. The first is the maximum level of refine- 
ment, which affects the spatial resolution and, ultimately, 
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Fig. 10. — Differential kinetic energy flux processed by shocks as a function of Mach number and pre-shock gas phase. Pre-shock gas 
overdensity (top row) is divided into several ranges that differentiate the overall distribution. Each pre-shock temperature range (bottom 
row) roughly corresponds to a particular mach number. Both distributions are shown for three redshifts of z = 3 (left column), z = 
1.5 (middle column) and 2 = (right column). 
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Fig. 11. — The effects of spatial resolution in different density 
regimes. Here we keep the mass resolution at the highest level 
(Mdm = 9-7 x 10 8 M Q ) . The distribution function of shock Mach 
numbers in a (32 Mpc/h) 3 volume around a cluster weighted by 
surface area for three overdensity cuts is plotted against Mach num- 
ber. Three cuts in overdensity are shown for 8 < 100 (black lines), 
100 < 5 < 10 4 (blue lines), <5 > 10 4 (red lines). Varying spatial res- 
olution are shown with dotted (62.4 kpc/h), dashed (15.6 kpc/h), 
and solid (3.9 kpc/h) linestyles. 

the accuracy of the hydrodynamic solver. The second 
parameter is the dark matter particle mass resolution, 
which affects the accuracy with which the gravitational 
potential is calculated. 

5.1. Spatial Resolution 

Our maximum spatial resolution ranges from 
62.4 kpc/h to 3.9 kpc/h (see Table [J). Since this 



Fig. 12. — Spatial resolution effects on kinetic energy flux for 
several density regimes. The kinetic energy flux is plotted against 
Mach number for varying spatial resolution denoted by linestyle. 
The shocks are grouped as external (black lines), clusters/filaments 
(blue lines), and rich clusters (red lines). 

only limits the maximum resolution, one expects to see 
a strong dependence on this parameter only at high 
densities. Figure [11] shows the dependence of shock 
surface area on level refinement for three overdensities. 
For shocks with pre-shock overdensities less than ~ 100, 
the main difference in the multiple resolutions is at low 
Mach numbers (below M ~ 2) and only appears in the 
lowest resolution simulation. 

At overdensities above 10 3 , we not only see that the 
low Mach number complex flow shocks are lost at low 
resolution, but also a drop in the number of high Mach 
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Fig. 13. — The effects of mass resolution in different density 
regimes. Here we keep the spatial resolution at the highest level 
of 3.9kpc/h. Mass resolutions are shown by linestyles of dotted 
(6.2 X 1O 1O M ), dashed (7.8 X 10 9 A/q), and solid (9.7 X 10 8 M Q ). 
The shocks are grouped as external (black lines), clusters/filaments 
(blue lines), and rich clusters (red lines). 

number shocks. At this density and spatial refinement, it 
is thought that the absence of sufficient spatial resolution 
leads to the artificial smoothing of the gas, creating an in- 
ability to capture shocks. The main result of this spatial 
resolution study is that a resolution between 3.9 kpc/h 
and 15.6 kpc/h should be sufficient in all but the most 
dense regions of the simulations. Therefore, our "Santa 
Fe Light Cone" simulation presented in Section 3 had an 
adequate peak resolution of 7.8 kpc/h. 

To examine the effect of spatial resolution on cosmic 
ray acceleration, we study the kinetic energy flux through 
shocks as a function of spatial resolution. Figure [T2] 
shows a distribution function measuring the thermal dis- 
sipation rate as a function of Mach number with vary- 
ing spatial resolution. This study is performed with the 
maximum mass resolution, M dm = 9.7 x 10 8 M o . At 
low overdensities (6 < 100), the effect of spatial reso- 
lution is very small. At moderate to high overdensities 
(100 < 5 < 10 4 ), there are differences on the order of a 
factor of 2 that are likely due to the smoothing of high 
density gas as the resolution is decreased. The primary 
difference in the dissipation rates occur for low Mach 
numbers when we do not have sufficient spatial resolu- 
tion to resolve all of the complex flow shocks. There are 
also large differences at M. > 10 for the lowest resolution 
simulation. However, the difference between 3.9 kpc/h 
and 15.6 kpc/h is negligible. 

At very high overdensities (5 > 10 4 ), there is a very 
large difference between the varying spatial resolutions. 
One reason is that if a cell has an overdensity of 10 , 
the grid would normally be on the 5th level of refine- 
ment. With a maximum refinement level of 4 for the 
poorest resolution simulation, any gas at this overdensity 
would be very poorly resolved. The difference between 
the 15.6 kpc/h and 3.9 kpc/h resolution simulations is 
likely small number statistics for the former simulation. 
The 3.9 kpc/h resolution simulation will resolve these 
high densities with roughly 64 times more cells compared 
to the 15.6 kpc/h simulation. 

5.2. Mass Resolution 

The mass resolution of each simulation is set by the 
resolution of the root grid (or highest-level static nested 
grid). The size of each root cell determines the amount 
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Fig. 14. — Spatial resolution effects on kinetic energy flux for 
varying dark matter particle masses. Mass resolutions are shown 
by linestyles of dotted (6.2 X 10 10 M Q ), dashed (7.8 X 1O 9 M ), and 
solid (9.7 X 10 s Mq). The shocks are grouped as external (black 
lines), clusters/filaments (blue lines), and rich clusters (red lines). 

of mass given to each dark matter particle. Therefore, if 
the root grid doubles in resolution, the mass resolution 
increases by a factor of 2 3 . In principle, there should be 
two effects of increased mass resolution. First, one would 
expect that since we are extending our mass function to a 
lower limit, the number of subhalos and our resolution of 
complex fluid flow should increase. This should manifest 
itself in an increase of shocks in the low Mach number 
regime. Second, the increased mass resolution also corre- 
sponds to an increase in the static grid spatial resolution. 
This may affect the calculated surface area of shocks that 
reside in voids. Since the temperature jumps in the voids 
are likely to be much higher than those within clusters, 
we would expect this effect to show up in the high Mach 
number regime. 

In order to test these hypotheses, we varied the mass 
resolution from 6.2 x 10 w M Q /h to 9.7 x 10 8 M Q /h. The 
results of this study are shown in Figure fl3l At 5 < 100, 
we see that as the mass resolution increases, the num- 
ber of low Mach number shocks increase, while the high 
Mach number shocks decrease. At high densities, the sit- 
uation is more complicated. For M. < 2, the surface area 
likely increases because of the increase in the number of 
subhalos and complex flow. For M > 7 — 8, the differ- 
ences seem to be largely due to statistical uncertainties. 
For 5 > 10 4 , the disparity at M < 2 is again likely due to 
the number of subhalos and their effects on turbulence. 
At 2 < M < 5, there is a large difference between the 
highest mass resolution simulation and the other two. 
Because we believe these shocks are merger shocks, it 
may be because there are just too few dense subhalos 
that merge with large halos to create this surface area. 

As with the spatial resolution, we now study the ef- 
fects of mass resolution on the thermal dissipation rates 
at shock fronts. Figure [14] shows the effect of var- 
ied mass resolutions with a fixed spatial resolution of 
3.9 kpc/h. Again, we break the analysis down into over- 
density regimes. Low overdensity, high Mach number 
shocks exhibit a strong dependence on the mass reso- 
lution. This is because of the ability to better resolve 
shocks in the voids and low density filaments. The dis- 
parity in high overdensity (100 < S < 10 4 ), low Mach 
number shocks is less apparent, but also suggests that 
the mass resolution of the simulation has an effect on the 
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thermal dissipation of gas through shocks. At 5 > 10 4 , 
we again see the effect of a decreased number of subhalos 
available to merge. 

Contrary to the effects of spatial resolution, the biggest 
differences due to mass resolution appear in the high 
Mach number regime. This is again due to the overesti- 
mate of Mach number at low root grid resolution. One 
evident result is that while the spatial resolution seems 
to be fairly well converged, it is not clear that the mass 
resolution has converged. Therefore, we can only claim 
a fairly weak precision in the thermal dissipation and 
cosmic ray acceleration rates for the current simulations. 

There are several key results to this resolution study. 
We appear to have converged in terms of maximum spa- 
tial resolution in all but the densest cluster gas. However, 
our convergence upon the various quantities with respect 
to dark matter mass resolution is not clear. The differ- 
ences in the cosmic ray acceleration rate are not larger 
than the underlying uncertainty in the results of diffu- 
sive shock acceleration simulations, suggesting that both 
mass resolution and our understanding of the physical 
mechanisms of cosmic ray acceleration must be improved 
in the future. 

6. DISCUSSION 

There are several topics that warrant discussion with 
respect to the results that we have presented thus far. 
These include the variation of results with respect to <7g, 
the inclusion of non-adiabatic physics, the limitation of 
the diffusive shock acceleration model, and the implica- 
tions of the mass resolution in the "Santa Fe Light Cone" 
Simulation. 

If our goal is to do large statistical studies of galaxy 
clusters, changing the value of as will have significant 
effects. First, a higher as will greatly increase the num- 
ber of massive clusters in a given volume. By comparing 
the ryul024 simulation with the "Santa Fe Light Cone," 
with values of a s of 0.8 and 0.9, respectively, we see 
that this increases the frequency and strength of the high 
Mach number shocks. Additionally, this should increase 
the amount of kinetic energy that is processed by shocks 
since mergers will be more frequent. 

In all of our simulations thus far we have only used adi- 
abatic phys i cs. Pr evious studies, such as those done by 
iKang et al.l (|2007f h have found that when including ra- 
diative cooling and star formation that the shock proper- 
ties are still governed primarily by gravitational physics 
and that additional physics have little effect on overall 
distr ibutions at scale s large r than ~ 100h~ 1 kpc. How- 
ever, |Pfrommer^Fal] (|2007f ) found that at smaller scales, 
on the inside of clusters, the cosmic ray contribution to 
the overall pressure is greatly increased with the inclusion 
of radiative cooling. Additionally, we currently adopt 
a temperature floor of 10 4 K because of the lack of an 
ionizing background. This should instead be done in a 
self-consistent manner. 

While we are using resu lts of recent d i ffusive shock 
acceleration simulations bv IKang fe Jones! ()2007l ). there 
are assumptions and limitations that may have an effect 
on our results. We assume that the magnetic field is 
parallel to the shock normal, which yields the largest ef- 
ficiency for accelerating cosmic rays. Any deviation from 
this will likely cause decreases in the overall efficiency of 
the shocks as particle accelerators. Additionally, for low 



Mach numbers, knowledge of the pre-shock composition 
is very important and can lead to orders of magnitude 
differences in the acceleration efficiency. Therefore, being 
able to track the cosmic ray pressure in "on-the-fly" cal- 
culations will allow us to provide a more self-consistent 
estimate. Finally, we are assuming that the only method 
for cosmic ray production is through first-order Fermi ac- 
celeration, and therefore we ignore other potetial sources 
of cosmic rays, such as second-order acceleration by tur- 
bulence, galaxies, and AGN. 

The results of the resolution study provided in Section 
4 have suggested that we have not yet seen a conver- 
gence with respect to the dark matter particle mass in 
the "Santa Fe Light Cone" simulation. This likely results 
in an under-prediction in the number of merging subha- 
los and the kinetic energy flux associated with them. 

7. CONCLUSIONS & FUTURE DIRECTIONS 

Our study of cosmological shocks has resulted in sev- 
eral advances in both scientific understanding and nu- 
merical algorithms. We now summarize the key findings: 

• We have developed a novel numerical scheme that 
is capable of detecting and accurately characteriz- 
ing the Mach number of shocks in an adaptive mesh 
refinement simulation. This method has relaxed 
the previous restriction of using a coordinated axis- 
based approach and now allows us to accurately 
characterize shocks that have any orientation with 
respect to the coordinate grid. 

• Using our new shock-finding technique on a uni- 
grid cosmological simulation tha t is identica l to th e 
highest-resolution calculation in lRyu et~a l. (2003), 
we have shown that previous methods resulted in 
an overestimate of the number low Mach number 
shocks by a factor of ~ 3 due to confusion of the 
direction of shock propagation, and that this un- 
derestimate is consistent with using shock-finding 
algorithms that only sweep along coordinate axes. 

• We have analyzed the largest AMR cosmologi- 
cal simulation to date that includes adiabatic gas 
physics, the "Santa Fe Light Cone." This simu- 
lation has an effective spatial dynamic range of 
65, 536 and resolves both large scale structure 
and small-scale features within galaxy clusters. 
Whereas previous studies were able to study on 
the order of 10 high mass clusters, we have thou- 
sands within a single simulation volume. Our study 
of this simulation has led to a new technique for 
conceptualizing structure formation because we are 
able to analyze the evolution of three different 
populations of shocks: cluster accretion, filament 
accretion, and internal merger and complex flow 
shocks. 

• By applying the results of 1-D Diffusive Shock Ac- 
celeration models, we calculate the amount of ki- 
netic energy at shock fronts that is used to accel- 
erate cosmic rays, and find it to be in agreement 
with previous studies. These cosmic rays will make 
up a significant fraction of the total pressure in the 
intracluster medium and therefore their dynamical 
effects need to be studied. 
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• We have performed a resolution study that varies 
both the dark matter particle mass and peak spa- 
tial resolution. From the results of this study, we 
believe that the spatial refinement in the "Santa 
Fe Light Cone" simulation is adequate. The state 
of mass resolution convergence is less clear, sug- 
gesting that for future studies we should focus on 
higher mass resolution. 

While our numerical technique of characterizing shocks 
has been proven to be robust, our results are still some- 
what limited by the physics. We have not yet included 
potentially important effects such as radiative cooling, 
star formation and feedback, AGN heating, or a pho- 
toionizing UV background. These physics will be in- 
cluded in future work. Second, our results are based 
upon a post-processing of the simulation output. Ideally, 
the shocks would be identified in an "on-the-fiy" manner 
during simulation runtime. Additionally, the cosmic ray 
acceleration would be traced in a self-consistent manner 
that allowed for a back-reaction on the gas. Attempts 
at tracing the c o smic r ay pr essure have been made by 
iPfrommer et all (|2006l . l2007f) using an SPH code, and 
we will be working towards the same goal in the near 
future within Enzo. Finally, the acceleration of cosmic 
rays is still dependent on the underlying magnetic field 



strength and orientation. Cosmological MHD has been 
implemented within Enzo, and in the near future we will 
include magnetic fields and their coupling to cosmic rays 
within a cosmological AMR volume. 
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